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PACS 05.45.Xt - Synchronization; coupled oscillators 

PACS 89 . 75 . Fb - Structures and organization in complex systems 

PACS 89 . 75 . He - Networks and genealogical trees 

Abstract. - In this Letter we investigate networks that have been optimized to realize a trade-off 
between enhanced synchronization and cost of wire to connect the nodes in space. Analyzing the 
evolved arrangement of nodes in space and their corresponding network topology a class of small 
world networks characterized by spatial and network modularity is found. More precisely, for low 
cost of wire optimal configurations are characterized by a division of nodes into two spatial groups 
with maximum distance from each other, whereas network modularity is low. For high cost of 
wire, the nodes organize into several distinct groups in space that correspond to network modules 
connected on a ring. In between, spatially and relationally modular small world networks are 
found. 
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Introduction. — Synchronization phenomena occur 
ih a diverse range of contexts in nature, engineering and 
society: cardiac pacemaker cells, neurons in the brain, 
fireflies that flash in unison, the power grid or consen- 
sus formation among people are just a few example ap- 
plications from these fields. All of these arc distributed 
systems embedded in space in which most couplings are 
local, but often also non-local long range couplings are 
present. Hence, most of these systems can be described as 
small world (SW) networks [I]: nodes represent elemen- 
tary units of the system such as neurons, fireflies, power 
stations or people and links in the network represent inter- 
actions that describe how the elementary units influence 
each other. Synchronization phenomena on SW and other 
complex networks have found much attention in the re- 
cent literature, see, e.g. [2], for a review. Moreover, in 
the system of coupled oscillators which we consider below, 
superior synchronization is essentially related to maximiz- 
ing the second smallest eigenvalue of the graph Laplacian. 
This eigenvalue -the algebraic connectivity- is an impor- 
tant invariant for undirected graphs and is, e.g., relevant 
for the analysis of the robustness of networks against node 
removal or for epidemic spreading [3J. 

In this research, one focus of interest has been to under- 
stand how the structure of the interaction network influ- 



ences the dynamics of synchronization. Various network 
properties have been connected with superior synchroniza- 
tion and, even though at least full synchronization ulti- 
mately depends on fine details of the network structure [5] , 
rough rules of thumb are that enhanced synchronization is 
positively correlated with short average pathlengths, dis- 
assortativity, large girths, and very homogeneous degree 
distributions. For instance, it has been found that syn- 
chrony optimal networks, termed 'entangled' networks by 
the authors of [4], are regular graphs. 

Among the research on synchronization characteristics 
of networks, the problem of optimal synchronization with 
constraints has found much interest. For instance, various 
studies reported about the weighted link arrangement that 
gives optimal synchronization on a given network topol- 
ogy 5 - 8 , the optimal connection architecture of strongly 
connected directed networks [5] or the structure of the 
'synchronization fitness landscape' in network space |10j . 
Other work addressed the problem of optimal synchroniza- 
tion of non-identical units on sparse networks |llU12j . 

As another constraint, the interplay between synchro- 
nization properties of a network and its embedding in 
space has recently been discussed [13] ■ In [T3J the spatial 
embedding of the network is represented as a constraint 
on the length of a wire needed to realize the connections of 
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the coupling network. It is found that in situations where 
the cost of these spatial connections is high, the network 
that gives optimal synchronization is a SW, in which the 
length of the connections in space are distributed accord- 
ing to a power law P(l) oc l~ a , if P(l) gives the probability 
that a randomly picked link connects oscillators of spatial 
distance The costlier the wire, the larger the exponent 
o. 

Modelling networked systems whose evolution is guided 
by considerations of function as well as spatial constraints 
has applications for a number of biological and technical 
systems. In the context of synchronization one may, for in- 
stance, think of neuronal networks for which synchroniza- 
tion plays an important role in the information processing. 
Further, neuronal networks are systems that have evolved 
and cost-considerations of long and short links have cer- 
tainly played a role in shaping their structure. Interest- 
ingly, a number of recent studies that attempt to entangle 
the large-scale organization of functional brain networks 
[21] have revealed a modular small- world organization - a 
network architecture that is very naturally explained by 
the model of spatially constrained synchrony-optimal net- 
works discussed below. 

In this Letter we extent the model of [13] by another 
degree of freedom. Different to previous work we do not 
consider nodes as having a fixed position in space, but al- 
low nodes to change locations. We do not only ask 'How 
many links docs a SW need to synchronize?', but also ask 
the question: 'What is the relative arrangement of nodes 
in space that allows for superior synchronization with lim- 
ited coupling?'. As we will discuss below, depending on 
constraints, this leads to yet another class of SWs which 
realize superior synchronization: relationally and spatially 
modular SWs. 

The Model. — In more detail, similar to [JHSJUI] we 

study identical synchronization in systems composed of N 
linearly coupled identical oscillators 

Si = f{si) + a {g(sj) - g(si)}. (1) 

3 

In the above equation = f(si) gives the dynamics of the 
individual oscillators, while Ay is the adjacency matrix 
of the coupling network, a the coupling strength and the 
function g describes the "inner" coupling of the oscillators. 
For simplicity we restrict the study to undirected net- 
works, i.e. symmetrical coupling matrices Ay . Depending 
on the details of the oscillators /, the coupling strength 
and architecture oAy and the coupling function g, in the 
system ([T]) identical synchronization can occur. Analysing 
the stability of the fully synchronized state /(s) =0 Pec- 
ora and Carroll [14] have derived a 'master stability func- 
tion' which relates the stability to the eigenvalues of the 
Laplacian matrix Gy = Sij An — A^ of the coupling 
network. Since we are interested in undirected graphs, 
the spectrum of G is real. Assuming that the networks 
are connected there is exactly one zero mode (which cor- 



responds to perturbations along the synchronization man- 
ifold) and the eigenvalues A^ , i = 1 , . . . , N of G may be 
labelled in ascending order = Ai < A2 < ... < A^. In 
Ref. [14] it has been shown that for a large class of oscil- 
lators the stability of the synchronized state is related to 
a small eigenratio e = A similar synchronization 

measure can be derived for the discrete analogue of Eq. 
(H|), synchronization on coupled map lattices, cf. [3J. 

The eigenratio is thus a measure for the stability of the 
fully synchronized state. Its independence of the details of 
the oscillator system allows it to define a measure for the 
synchronization properties of a network. For this reason 
it has been widely adopted in various studies [5 HTU1I13II14] , 
there being also some support that it even gives a rough 
measure for synchronization properties of systems of non- 
identical oscillators [TJE] - 

To proceed, we introduce the spatial aspect of the sys- 
tem. We consider oscillators with locations Axi, i = 
1,...,N which are distributed on a 1-dimensional space 
with periodic boundary conditions. The parameter A de- 
fines the length of a measurement unit on that space, 
which we set to A = 1 in the following considerations. 
Defining D 

max — mflX^ Xi cl distance metric may then be 
defined by d(i,j) = min(|a;j — Xj\, D max — \%i — Xj\), such 
that the amount of wire needed to connect the oscillators 
is given by W = J2 t<j A tJ d(i,j). 

In the following, we are interested in network config- 
urations and spatial node arrangements that allow for 
superior synchronization properties while minimizing the 
amount of wire needed to realize them, i.e. networks that 
minimize an 'energy' 

E = (3W/N + (1 - (3)e. (2) 

In Eq. (J2J) the parameter j3 determines the relative de- 
sirability of both factors, superior synchronization prop- 
erties (measured by the eigenratio e) and minimal cost of 
connections (measured by the amount of wire per node 
W/N). The formalism of investigating trade-offs in net- 
work formation in Eq. (0) is very similar to approaches 
in previous studies, e.g., [13j[T5l[T6]. As it stands, the 
solution to the problem @ is trivial if the locations of 
the oscillators are not fixed in space: all oscillators would 
move to one location, thus allowing for full connectivity 
at no cost of wire. In a real-world situation constraints 
such as limited mobility of the nodes, physical barriers to 
movement, minimum space requirements of nodes or other 
functional reasons that limit nodes to certain parts of the 
space prevent this configuration. We model the sum of 
these 'real-world' limitations as a constraint that the av- 
erage spatial distance D(x) = X^j<j between nodes 
is held constant during the minimization of ([2]). The for- 
mulation of the problem in this framework allows to study 
what relative arrangement of oscillators in space gives rise 
to optimal synchronization for least cost of wire. 

We continue by numerically constructing oscillator con- 
figurations that minimize E(f3) subject to D(x) = const. 
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Fig. 1: Examples of evolved networks for different trade-offs between cost of wire and desireability for superior synchronization: 
(a) p = O.Of (very low cost of wire), (b) /3 = 0.5 (balanced costs for wire and synchronization), and (c) /3 = O.Of (very high cost 
of wire). The networks are of size N = fOO and contain L — 400 links. In the figure vertices have been colored according to the 
modules they belong to (modularities are Q = .26 for (a) and Q — 0.71 and Q = 0.78 for (b) and (c)). The spatial locations 
roughly correspond to the evolved spatial locations of the nodes during the optimization, however a random number was added 
to make vertices distinguishable. 



for different trade-off parameters (3. The numerical opti- 
mization scheme, which implements a variant of simulated 
annealing, is similar to schemes employed in previous stud- 
ies, like, e.g., (U[T3J[T5J[TB]. It essentially consists of the 
following steps: 

(i) Start with a configuration of oscillators that arc 
evenly distributed on a 1-dimensional ring, i.e. x® = i, 
i = 1, N, which also defines the average spatial distance 
of the oscillators Dq = D(x ). The coupling matrix in the 
initial condition is assumed to be given by an Erdds-Renyi 
random graph [T7] with exactly L links, (ii) Suggest either 
a rewiring or a location change for one or several oscilla- 
tors. For rewiring suggestions, moves of randomly selected 
links to randomly selected 'link vacancies' are suggested 
(provided that the new links do not introduce self-loops 
and that the new configuration is connected) . In the case 
of a location change, a move of a randomly selected oscil- 
lator location Xi — ¥ Xi + Sxi (provided that Xi + Sxi > 0) 
and rescaling Xi — > XiD(x) / D(x + Sx) of the oscillator lo- 
cations such that D(x) = D is suggested, (iii) Calculate 
the energy E'(/3) of the modified configuration and accept 
if E' < E or with probability p = exp(— v(E' — E)) oth- 
erwise. As usual in simulated annealing procedures the 
'inverse-temperature parameter' v is gradually increased 
during the optimization. If, according to the above rule, a 
suggested configuration is not accepted, the previous con- 
figuration is restored, (vi) Iterate steps 2 and 3 until no 
improvement in E{f3) could have been obtained for the 
last 10L iteration steps. The motivation for the terminat- 
ing condition is to stop the algorithm after every link or 
node location has been unsuccessfully tried to be modified 
several times. 

It should be emphasized that the above optimization 
problem is a difficult non-linear problem and the numeri- 
cal procedure does not ensure that a global optimum has 
been reached. However, we repeated the stochastic ex- 



periment for different initial conditions and the features 
about optimal configurations that we report below have 
been found to be robust. 

Results. — Before discussing how varying trade-offs 
between synchronizability and cost of wire influence the 
structure of the optimal network configurations, it is 
worthwhile to investigate the limiting cases (3 = and 
/3 = 1 of our model. The case of (3 = 0, no cost of wire, 
i.e. when the spatial arrangement of nodes becomes irrele- 
vant, corresponds to the model of [4], which we briefly dis- 
cussed in the introduction. Importantly, optimal networks 
for (3 = are not modular and -since space is irrelevant - 
the spatial arrangement of nodes in the optimal configu- 
ration is uniform. The optimal configuration for the other 
limiting case f3 = I, minimization of the cost of wire with- 
out regard for synchronization, is an arrangement of nodes 
into two distinct spatial clusters of nodes at maximum dis- 
tance. Nodes within each of the two clusters are strongly 
interconnected, such that they also correspond to two net- 
work modules. These two modules are connected by ex- 
actly one link. Due to its community organization this is 
a network configuration with very poor synchronizability, 
cf. Ql]. 

To proceed, Figure Q] shows some representative exam- 
ple networks of N = 100 nodes and L = 400 links evolved 
for three different values of the trade-off parameter (3: an 
example when wire is inexpensive for (3 = 0.01 (cf. Fig. 
QJ,), an example for when the cost of wire and desirabil- 
ity of superior synchronization properties are balanced for 
(3 = 0.5 (cf. Fig. [T})) and an example for a situation when 
wire is very expensive for f3 = 0.99 (cf. Fig. \Tjp). In the 
figures, nodes with spatially close locations are also drawn 
close to each other, even though distances are sometimes 
magnified or shrunk to make individual nodes distinguish- 
able. The illustrations of the synchrony-optimized exam- 
ple networks strongly suggest the emergence of different 
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Fig. 2: Dependence of some properties of the evolved networks 
and spatial arrangements on the trade-off parameter /3: (a) 
eigenratio e = (b) average link length (measured in 

units of the maximum coordinate), (c) clustering coefficient C, 
(d) modularity Q. For reference, the horizontal lines indicate 
the range the respective quantities would assume for an Erdos- 
Renyi random graph whose nodes are uniformly distributed in 
space. In (b) the lines are omitted for scaling reasons, one has 
W/(LD max ) w 0.5. 



Fig. 3: Dependence of (a) the average shortest pathlength d, 
(b) average diameter d ma x, (c) spatial modularity S, and (d) 
average relative spatial distance between nodes in the same 
module in units of D ma x on the trade-off parameter /3. For 
reference, the horizontal lines indicate the range the respec- 
tive quantities would assume for an Erdos-Renyi random graph 
whose nodes are uniformly distributed in space. In (c) and (d) 
the lines are omitted for scaling reasons, one has S < 1CF 3 and 
^mod ~ -25, respectively. 



types of modular network organizations that depend on 
the trade-off parameter (3. 

To explore the relational or network modularity, we have 
calculated the modularity Q = ^2 m [L m /L — (d m /2L) 2 } 
introduced in [19 for the networks. In the expression for 
Q the index m extends over all modules, L m denotes the 
number of links between nodes of module to, and d rn is the 
sum of the degrees of all nodes in module to. Since the 
example networks investigated are relatively small we have 
analyzed the modules via extremal stochastic optimization 
[20) . but have also tried other methods like [19], which 
have robustly confirmed the visual expectation one gathers 
from the network plots in Fig. [TJ in which nodes belonging 
to the same network module are characterised by identical 
colouring. 

Clearly, when wire is inexpensive and superior synchro- 
nization the dominant consideration, nodes assemble into 
two spatial clusters that are separated by the maximum 
spatial distance. In this case, however, the network- or re- 
lational structure of the coupling has very low modularity 
and the modules found by the algorithm appear uncorre- 
cted with the spatial arrangement to the eye. 

When the demands for wire minimization and superior 
synchronization are balanced networks like the one visu- 
alized in Fig. [TJ} emerge. The example network clearly 



decays into several distinct relational modules and al- 
ready has a high network modularity of Q = 0.71. Like- 
wise, however, spatial clusters of nodes emerge in close 
correspondence with the network modules. Even though 
the modules are clearly distinct, they are still relatively 
strongly interconnected, mostly by long distance links. 

The tendency towards stronger spatial and relational 
modularity continues when (3 is further increased. For f3 = 
0.99, when wire economy is the main consideration in the 
optimization, very cohesive spatial and network modules 
can be discerned, cf. Fig. [JJ:. Unlike as for the balanced 
case, these modules are connected to each other in a ring 
like arrangement by short links between spatial nearest 
neighbour modules. 

For a more systematic investigation we constructed 
R = 100 optimized network configurations for system- 
atically varied trade-off parameters f3. All networks in- 
vestigated below are of size N = 400 and have L = 400 
links. By displaying some key statistics of network ar- 
rangement and node arrangement in space, figures [2] and 
[3] give an overview over the parameter space. By plotting 
the eigenratio A2 / \n and the average link length (in units 
of the maximum spatial coordinate Z? max ), panels (a) and 
(b) of Fig. [5] visualize the trade-off between superior syn- 
chronization properties and cost of wire. As expected, 
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for small /3 much wire can be used to achieve superior 
synchronizabilities, whereas for large (3 costly wire allows 
for only few long links, poor synchronization behaviour 
being the consequence. To further classify the optimal 
network topologies, we investigated their degree distribu- 
tions, degree variance, clustering coefficients (as defined 
in [I]), modularity, and the average shortest pathlengths 
and diameters. In agreement with we find that for 
the full range of (3 parameters the degree distributions of 
the evolved networks are very narrow and the degree vari- 
ances are significantly lower than for random graphs with 
the same number of nodes and links. However, as reported 
in |I3j , degree variances do not become exactly zero when 
/3 > and increase when f3 is increased, i.e. when spatial 
constraints become more important. 

More interestingly, panels (c) and (d) of Fig. [^illustrate 
that the optimized networks become increasingly cliquish 
and modular when the relative cost of wire is increased. 
Thus, even for balanced cost of wire and desirability of 
superior synchronization, when synchronizablities in the 
range of Erdos-Rcnyi random graphs are attained, the net- 
work organization is characterized by a distinctly modular 
arrangement, cf. also Fig. [Ip. While for low (3 the mod- 
ules are connected by many long range links, these long 
range connections are increasingly thinned out when j3 in- 
creases. As a consequence, while the modular arrangement 
becomes more and more distinct, diameter and shortest 
pathlengths of the optimized networks grow, cf. panels 
(a) and (b) of Fig. [3] 

On the other hand, when wire is cheap, triangles are 
suppressed and no relational community organization with 
more cohesiveness than found in a random graph is de- 
tected, cf. Fig. [5£,d. In this case the networks also be- 
come increasingly smaller, approaching the entangled net 
configurations discussed in [3]. 

As the network plots in Fig. [T] already show, the opti- 
mized configurations are not only characterized by a dis- 
tinct network arrangement, but also have a characteristic 
arrangement of nodes in space. To classify the spatial 
arrangement of the nodes, we calculated a spatial correla- 
tion function G(x) that measures the average probability 
of finding a node at spatial distance x from an arbitrary 
node. Figure H] displays plots of G(x) for two typical spa- 
tial node arrangements found for (3 = 0.05 and (3 = 0.99. 
The presence of one or two distinct peaks in the function 
indicates a high degree of spatial clustering. While the two 
peaks comprising each roughly N/2 nodes for (3 = 0.05 
confirm the presence of two spatial clusters separated by 
the maximum spatial distance, the one peak that com- 
prises about N/6 nodes for (3 = 0.99 indicates the presence 
of several spatial clusters. In the multiple cluster config- 
uration, individual clusters are not separated by a clearly 
defined typical spatial length scale, but tend to avoid each 
other (i.e. G{x) has a broad minimum at ranges x = 0.02 
to a: w 0.12). The typical width of the peaks Ax = 0.01 
allows to define a more condensed spatial modularity mea- 
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sure S = J Q x G(x)dx, that quantifies the average number 
of nodes in the immediate spatial vicinity of an average 
node, thus defining the average relative size of a spatial 
module. Panel (c) of Fig. Ogives the dependence of 5" on 
(3 for the optimized network configurations. 

It is also of interest, whether spatial clusters overlap 
with network modules. To characterize this relationship, 
we calculated an overlap parameter Z mo d that measures 
the average spatial distance between nodes belonging to 
the same network module in units of the maximum co- 
ordinate Umax- If spatial and relational modules overlap 
completely, one expects to have / mo d < Aa;/D max , whereas 
'mod ~ 1/4 if spatial and relational modules are completely 
uncorrelatcd. The dependence of 'mod on (3 is plotted in 
panel (d) of Fig. H 

As discussed above, for j3 = 1 a network configuration 
characterized by two spatial clusters corresponding to two 
network modules is optimal. In this case one has S = 1/2 
and / m od ~ 0. Introducing a small consideration of syn- 
chronizability into the fitness function defined in Eq. ^ 
leads to a sharp drop in S, the network is immediately dis- 
tributed into a number of smaller modules that comprise 
many less than N/2 nodes. In this process the overlap 
of spatial and relational module structures is preserved. 
This almost complete overlap persists over a wide range 
of (3 parameters until around (3 ~ 0.07. This value of j3 
marks a sharp transition in the structure of typical opti- 
mized configurations. Below (3 = 0.07 two spatial clusters 
that do not overlap with network cliques are found. In 
this situation the link length distributions are marked by 
two sharp peaks: very short links connecting nodes in the 
same spatial cluster and many links of maximum length 
connecting nodes pertaining to different spatial clusters. 

For (3 = 0, the cost of wire and spatial constraints be- 
come irrelevant and thus nodes in optimal network config- 



p-5 



M. Brcdc 



urations are distributed uniformly in space. Accordingly, 
below the transition at (3 = 0.07 the spatial modularity 
first sharply increases when the spatial two cluster con- 
figuration becomes established. However, after reaching 
a maximum at ft = 0.05 the spatial modularity declines 
again, approaching S = Ax/Z? max from above. 

Discussion and conclusions. In this Letter we 
have numerically constructed spatially embedded optimal 
networks that realize a trade-off between superior synchro- 
nization properties and cost of wire. Different from all pre- 
vious work on optimal synchronization, we focussed on the 
interplay between the relational- and spatial organization 
of the optimized networks. 

As the most interesting feature of our analysis we find 
that the optimal configurations arc characterized by an in- 
terplay of spatial clustering and network modularity. Es- 
sentially two parameter regimes, which are separated by 
a sharp transition, have been identified. When considera- 
tions of enhanced synchronizability outweigh requirements 
for the economy of wire (0 < j3 < 0.07), the optimal con- 
figurations are characterized by an arrangement of nodes 
into two spatial clusters with many long range connections 
between them, but also a remainder of short connections 
which link nodes within the same spatial cluster. How- 
ever, the link arrangement is such that the spatial clusters 
do not correspond to relational (or network-) modules. 
When wire economy is important but not dominant for 
0.07 < p < 1, the optimal networks were found to de- 
cay into a number of clearly separated network modules. 
These network modules closely correspond to spatial clus- 
ters of nodes. 

Our work has a number of interesting implications. 
First, it suggests a new explanation for the emergence of 
network modularity via the minimization of the cost of a 
wire when nodes are free to arrange themselves in space 
during the minimization procedure. A spatial arrange- 
ment of nodes in clusters serves to reduce the cost of wire: 
links between nodes in the same spatial cluster can essen- 
tially be introduced without or with very little cost. The 
interplay of this process with an additional mechanism 
which favours an unmodular network arrangement (in our 
case the demand for enhanced synchronization) can cause 
the breakup of large modules and result in a network ar- 
rangement characterized by the presence of a large number 
of small network modules. There is a variety of candidates 
for such processes, the simplest of which is probably the 
minimization of average shortest pathlcngths, which has 
been considered in [T5J[TB]. 

Second, in the context of synchronization processes on 
networks, our work demonstrates that an additional con- 
straint -minimization of the cost of a wire needed to con- 
nect nodes embedded in space- which appears plausible 
in the context of biological or technical applications can 
cause network organizations to be 'optimal', which have 
before been identified as suboptimal when synchronization 
properties are analyzed by purely investigating the struc- 



ture of the coupling network. Dynamically, this modu- 
lar structure is associated with the presence of different 
timcscales for synchronization |18j . Our work suggests 
that oscillators which synchronize at the same timescale 
will also be located close to each other in space, a finding 
that may have applications in neurobiology [21j . 

Third, it appears of interest that when considering 
nodes that are free to move in space, power laws in the link 
length distribution of optimal networks reported in [13, 1G 
are replaced by bimodal link length distributions. Hence, 
an investigation of link length distributions of, e.g., bio- 
logical networks may allow conclusions about the nature 
of evolutionary processes that shape the networks. 

This research was undertaken on the NCI National Fa- 
cility in Canberra, Australia, which is supported by the 
Australian Commonwealth Government. I thank A. Kallo- 
niatis and F. Boschctti for helpful comments. 
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